### Configuration
################################################################################
### R Options
nodeInfo = system(paste0("scontrol show node ", Sys.info()["nodename"]), intern = TRUE)
cpus = as.integer(unlist(strsplit(unlist(strsplit(nodeInfo[grep("CfgTRES", nodeInfo)], "cpu="))[2], ","))[1])
memInKB = as.integer(unlist(strsplit(unlist(strsplit(nodeInfo[grep("CfgTRES", nodeInfo)], "mem="))[2], "M,"))[1]) * 1024^2
options(stringsAsFactors=FALSE,
dplyr.summarise.inform=FALSE,
future.globals.maxSize=min(memInKB, 20 * 1024^3),
mc.cores=min(cpus,1),
future.fork.enable=TRUE, future.plan="multicore",
future.rng.onMisuse="ignore")
### Fuctions
# Git directory and files to source must be done first, then all helper functions can be sourced
git_files_to_source = c("functions_analysis.R",
"functions_biomart.R",
"functions_degs.R",
"functions_io.R",
"functions_plotting.R",
"functions_util.R")
git_files_to_source = file.path(param$path_to_git, "R", git_files_to_source)
file_exists = purrr::map_lgl(git_files_to_source, file.exists)
if (any(!file_exists)) stop(paste("The following files could not be found:", paste(git_files_to_source[!file_exists], collapse=", "), ". Please check the git directory at '", param$path_to_git, "'.!"))
invisible(purrr::map(git_files_to_source, source))
### Debugging mode
# 'default_debugging' for default, 'terminal_debugger' for debugging without X11, 'print_traceback' for non-interactive sessions
param$debugging_mode = "default_debugging"
switch (param$debugging_mode,
default_debugging=on_error_default_debugging(),
terminal_debugger=on_error_start_terminal_debugger(),
print_traceback=on_error_just_print_traceback(),
on_error_default_debugging())### Rmarkdown configuration
################################################################################
### Create output directories
if (!file.exists(file.path(param$path_out, "figures"))) dir.create(file.path(param$path_out, "figures"), recursive=TRUE, showWarnings=FALSE)
if (!file.exists(file.path(param$path_out, "data"))) dir.create(file.path(param$path_out, "data"), recursive=TRUE, showWarnings=FALSE)
### R Options
options(citation_format="pandoc",
knitr.table.format="html",
kableExtra_view_html=TRUE)
### Knitr default options
knitr::opts_chunk$set(echo=TRUE, # output code
cache=FALSE, # do not cache results
message=TRUE, # show messages
warning=TRUE, # show warnings
results = "hold", # show results after running whole chunk
tidy=FALSE, # do not auto-tidy-up code
fig.width=10, # default fig width in inches
class.source='fold-hide', # by default collapse code blocks
dev=c('png', 'svg', 'tiff'), # create figures in png, tiff, and svg; the first device (png) will be used for HTML output
dev.args=list(png = list(type="cairo"), # png: use cairo - works on cluster
tiff = list(type="cairo", compression = 'zip')), # tiff: use cairo - works on cluster, transparent background, compression type "zip" is ‘deflate (Adobe-style)’
dpi=300, # figure resolution
fig.retina=2, # retina multiplier
fig.path=paste0(file.path(param$path_out, "figures"), "/") # Path for figures in png and pdf format (trailing "/" is needed)
)
### Set output hooks
knitr::knit_hooks$set(message=format_message, warning=format_warning)
### Required libraries
library(magrittr)
### Set the bibliographic environment
# Clear previous bibliography
knitcitations::cleanbib()
# If the DOI publication servers cannot be reached, there will be no citations, knitcitations will not write a references.bib file and pandoc will stop. This makes sure that there is always at least one citation.
bib_ref = knitcitations::read.bibtex(file.path(param$path_to_git,"assets","sc_analysis_references.bib"))
invisible(knitcitations::citep(bib_ref))
### Figure heights
# high figures
fig_large_height = 8
# Single figure, landscape
fig_standard_height = 4
# Two plots alongside (e.g. umaps)
fig_standard2_height = 5
# Three plots alongside (e.g. umaps)
fig_standard3_height = 4
# Four plots alongside (e.g. umaps)
fig_standard4_height = 2.5
# Four plots 2x2 (e.g. umaps)
fig_patchwork4_height = fig_standard2_height * 2
# Six plots 2x3 (e.g. umaps)
fig_patchwork6_height = fig_standard3_height * 2
# Eight plots 4x2 (e.g. umaps)
fig_patchwork8_height = fig_standard4_height * 2
# Twelve plots 4x3 (e.g. umaps)
fig_patchwork12_height = fig_standard4_height * 3
# Sixteen plots 4x4 (e.g. umaps)
fig_patchwork16_height = fig_standard4_height * 4# Load renv
renv::load(file.path(param$path_to_git,"env/basic"))
# Required libraries
library(liana)
library(ComplexHeatmap)
library(patchwork)
library(magrittr)
library(ggplot2)
library(ggsci)Cell-cell communication (CCC) is a process by which cells react to stimuli during many biological processes. CCC commonly focuses on ligand binding to its corresponding receptors but also includes interactions via extra-cellular matrix proteins and transporters. The stimuli commonly elicit in the receiver cells (the cells receiving the signal) a downstream response such as the induction of canonical pathways and transcription factor binding.
CCC inference methods usually use one of the following approaches: 1) Prediction of CCC interactions alone, commonly referred to as ligand-receptor inference as performed by LIANA 2) Analysis of downstream intracellular activities induced by CCC as conducted by NicheNet
Here, we will use LIANA to infer the ligand-receptor interactions as described in the LIANA tutorial (https://saezlab.github.io/liana/articles/liana_tutorial.html) and the ‘LIANA x Tensor-cell2cell Quickstart’ vignette (https://ccc-protocols.readthedocs.io/en/latest/notebooks/ccc_R/QuickStart.html).
The results of all CCC inference methods are affected by ligand-receptor prior knowledge resources i.e. the used databases (quality and contend, and each resource comes with its own biases). In this regard, LIANA provides a consensus expert-curated resource that combiness different resources (CellPhoneDBv223, CellChat27, ICELLNET30, connectomeDB202025, CellTalkDB31, as well as 10 others).
LIANA performs multiple CCC inference methods using the consensus resource and combines the results. The different CCC scores obtained by the diverse methods are categorized and translated into a “Magnitude” and “Specificity” value for the specific interactions. The “Magnitude” of an interaction is a measure of the strength of the expression of the interaction partners, while the “Specificity” indicates how specific an interaction is to a given pair sender and reciever cell tpye. Generally, both values are complementary. For example, a ligand-receptor interaction with a high magnitude score in a given pair of cell types is likely to also be specific.
We load the dataset as S4 class Seurat objects.
### Read rds object
################################################################################
### Load Seurat S4 objects
# Test if file is defined
if (is.null(param$data)) {
message("Dataset is not specified")
} else {
# Test if file exists
if (file.exists(param$data)) {
# Read object
message(paste0("Load dataset:", param$data))
sc = base::readRDS(param$data)
# Transfer original params to loaded object
if ("parameters" %in% list_names(sc@misc[])) {
# Retrieve previous parameter settings
orig_param = sc@misc$parameters
if ("SCT" %in% names(sc@assays)) {
if ("scale.data" %in% Layers(sc[["SCT"]])) {
orig_param$norm = "SCT"
} else {
orig_param$norm = "RNA"
}
} else {
orig_param$norm = "RNA"
}
# Keep some parameter settings from object and project defined
orig_param_keep = orig_param[c("annot_version", "species")]
basic_param_keep = param[c("path_to_git", "scriptname", "author", "project_id", "data", "path_out", "file_annot", "file_cc_genes")]
# Integrate parameter
param = modifyList(x = param, val = orig_param)
param = modifyList(x = param, val = basic_param_keep)
param = modifyList(x = param, val = param_advset)
param = modifyList(x = param, val = orig_param_keep)
}
### Set colors
# Set sample colors based on orig.ident
if (!is.null(sc@meta.data[["orig.ident"]])) {
if (length(unique(sc@meta.data[["orig.ident"]]))!=length(param$col_samples)) {
message(paste0("No or wrong number of distinct colors for samples provieded. Generate colors using color palette", param$col_palette_samples))
sample_colours_out = suppressWarnings(SetSampleColours(sc, param$col_palette_samples))
sc = sample_colours_out[[1]]
param$col_samples = sample_colours_out[[2]]
}
}
# Set colors for clusters based on seurat_clusters
if (!is.null(sc@meta.data[["seurat_clusters"]])) {
if (length(unique(sc@meta.data[["seurat_clusters"]]))!=length(param$col_clusters)) {
message(paste0("No or wrong number of distinct colors for clusters provieded. Generate colors using color palette", param$col_palette_clusters))
cluster_colours_out = suppressWarnings(SetClusterColours(sc, param$col_palette_clusters))
sc = cluster_colours_out[[1]]
param$col_clusters = cluster_colours_out[[2]]
}
}
# Set colors for cell types based on annotation
if (!is.null(sc@meta.data[["annotation"]])) {
if (length(unique(sc@meta.data[["annotation"]]))!=length(param$col_annotation)) {
message(paste0("No or wrong number of distinct colors for cell tpye annotation provieded. Generate colors using color palette", param$col_palette_annotation))
annotation_colours_out = suppressWarnings(SetAnnotationColours(sc, param$col_palette_annotation))
sc = annotation_colours_out[[1]]
param$col_annotation = annotation_colours_out[[2]]
}
}
sc
} else {
message("Dataset does not exist")
}
}×
(Message)
Load
dataset:/mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/output/Testdata/cluster_analysis/data/sc.rds
### Load reference Seurat S4 objects if specified
# Test if file is defined
if (!is.null(param$refdata)) {
# Test if file exists
if (file.exists(param$refdata)) {
# Read object
message(paste0("Load dataset:", param$refdata))
scR = base::readRDS(param$refdata)
# Transfer original params to loaded object
if ("parameters" %in% list_names(scR@misc[])) {
orig_paramR = scR@misc$parameters
if (!is.null(orig_paramR$col_samples)) {
param$col_samples_ref = orig_paramR$col_samples
}
if (!is.null(orig_paramR$col_clusters)) {
param$col_clusters_ref = orig_paramR$col_clusters
}
if (!is.null(orig_paramR$col_annotation)) {
param$col_annotation_ref = orig_paramR$col_annotation
}
param = modifyList(x = param, val = param_advset)
}
### Set colors
# Set sample colors based on orig.ident
if (!is.null(scR@meta.data[["orig.ident"]])) {
if (length(unique(scR@meta.data[["orig.ident"]]))!=length(param$col_samples_ref)) {
message(paste0("No or wrong number of distinct colors for samples provieded. Generate colors using color palette", param$col_palette_samples))
sample_colours_out = suppressWarnings(SetSampleColours(scR, param$col_palette_samples))
scR = sample_colours_out[[1]]
param$col_samples_ref = sample_colours_out[[2]]
}
}
# Set colors for clusters based on seurat_clusters
if (!is.null(scR@meta.data[["seurat_clusters"]])) {
if (length(unique(scR@meta.data[["seurat_clusters"]]))!=length(param$col_clusters_ref)) {
message(paste0("No or wrong number of distinct colors for clusters provieded. Generate colors using color palette", param$col_palette_clusters))
cluster_colours_out = suppressWarnings(SetClusterColours(scR, param$col_palette_clusters))
scR = cluster_colours_out[[1]]
param$col_clusters_ref = cluster_colours_out[[2]]
}
}
# Set colors for cell types based on annotation
if (!is.null(scR@meta.data[["annotation"]])) {
if (length(unique(scR@meta.data[["annotation"]]))!=length(param$col_annotation_ref)) {
message(paste0("No or wrong number of distinct colors for cell tpye annotation provieded. Generate colors using color palette", param$col_palette_annotation))
annotation_colours_out = suppressWarnings(SetAnnotationColours(scR, param$col_palette_annotation))
scR = annotation_colours_out[[1]]
param$col_annotation_ref = annotation_colours_out[[2]]
}
}
scR
} else {
message("Reference dataset does not exist")
}
} ## An object of class Seurat
## 6439 features across 1078 samples within 1 assay
## Active assay: RNA (6439 features, 3000 variable features)
## 3 layers present: data, counts, scale.data
## 2 dimensional reductions calculated: pca, umap
LIANA calculates an aggregate_rank across the used methods and generates a probability distribution for ligand-receptors that are ranked consistently better than expected under a null hypothesis. By default, the aggregate_rank in ascending order and can be interpreted similar to a p-value. Alongside we provide the aggregated “Magnitude” (named magnitude_rank) and “Specificity” (specificity_rank) of each interaction.
# Change omnipathR log level (ERROR or OFF)
OmnipathR::omnipath_set_loglevel(logger::ERROR, target = 'console')
liana_res = liana_wrap(sce = sc, idents_col = "annotation", assay = "RNA", assay.type = 'logcounts', method = param$liana_methods, resource = c('Consensus'), expr_prop = 0.1, permutation = list(nperms = 500), return_all = FALSE, verbose = FALSE)
# aggregate_rank
liana_agg = liana_aggregate(liana_res, verbose = FALSE)
write.csv(liana_agg, file.path(param$path_out, "data", "liana_aggregate_rank.csv"))
# magnitude_rank and specificity_rank
liana_rank_agg = rank_aggregate(liana_res, verbose = FALSE)
write.csv(liana_rank_agg, file.path(param$path_out, "data", "liana_results.csv"))
# aggregate_rank
knitr::kable(head(liana_agg, 25), align="l", caption="LIANA aggregate_rank, top 25 rows") %>%
kableExtra::kable_styling(bootstrap_options=c("striped", "hover")) %>%
kableExtra::scroll_box(width="100%")
# magnitude_rank and specificity_rank
knitr::kable(head(liana_rank_agg, 25), align="l", caption="LIANA magnitude_rank and specificity_rank, top 25 rows") %>%
kableExtra::kable_styling(bootstrap_options=c("striped", "hover")) %>%
kableExtra::scroll_box(width="100%")| source | target | ligand.complex | receptor.complex | aggregate_rank | mean_rank | logfc.logfc_comb | logfc.rank | natmi.edge_specificity | natmi.rank | cellphonedb.pvalue | cellphonedb.rank |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Monocytes | B-cells | TNFSF13B | TNFRSF13C | 0.0012215 | 47.00000 | 2.7856216 | 6 | 0.8761565 | 3 | 0 | 132 |
| Monocytes | B-cells | TNFSF13B | CD40 | 0.0027377 | 48.33333 | 2.3482866 | 9 | 0.8179126 | 4 | 0 | 132 |
| Monocytes | Monocytes | GRN | SORT1 | 0.0033754 | 49.00000 | 2.1961055 | 10 | 0.7359849 | 5 | 0 | 132 |
| Monocytes | Monocytes | PSAP | LRP1 | 0.0056819 | 50.66667 | 2.7616439 | 7 | 0.5199923 | 13 | 0 | 132 |
| Monocytes | Monocytes | PSAP | SORT1 | 0.0065810 | 51.33333 | 2.3709823 | 8 | 0.5100402 | 14 | 0 | 132 |
| Monocytes | B-cells | TNFSF13B | TNFRSF13B | 0.0065810 | 49.33333 | 1.9835360 | 14 | 0.9117102 | 2 | 0 | 132 |
| Monocytes | Monocytes | CD14 | TLR6 | 0.0108211 | 55.66667 | 1.9471684 | 17 | 0.4567520 | 18 | 0 | 132 |
| Monocytes | NK cells | IL1B | ADRB2 | 0.0108211 | 52.33333 | 1.9252123 | 18 | 0.5757635 | 7 | 0 | 132 |
| Monocytes | Monocytes | GRN | TNFRSF1B | 0.0120408 | 52.00000 | 2.8381942 | 5 | 0.4409884 | 19 | 0 | 132 |
| B-cells | CD8+ T-cells | HLA-DQA1 | LAG3 | 0.0146700 | 54.00000 | 1.8050511 | 21 | 0.5356698 | 9 | 0 | 132 |
| Monocytes | CD8+ T-cells | S100A8 | CD69 | 0.0175097 | 62.66667 | 3.7367472 | 1 | 0.2472976 | 55 | 0 | 132 |
| Monocytes | Monocytes | NRG1 | MS4A4A | 0.0175097 | 81.33333 | 0.6792135 | 111 | 0.9766946 | 1 | 0 | 132 |
| Monocytes | Monocytes | LGALS9 | LRP1 | 0.0206797 | 57.66667 | 1.6619519 | 25 | 0.4775093 | 16 | 0 | 132 |
| B-cells | CD8+ T-cells | HLA-DQB1 | LAG3 | 0.0223371 | 58.33333 | 1.6331241 | 26 | 0.4735673 | 17 | 0 | 132 |
| Monocytes | NK cells | IL1B | SIGIRR | 0.0240560 | 60.66667 | 1.7099145 | 23 | 0.3574421 | 27 | 0 | 132 |
| Monocytes | NK cells | SECTM1 | CD7 | 0.0240560 | 58.00000 | 1.5863438 | 27 | 0.4882453 | 15 | 0 | 132 |
| Monocytes | Monocytes | ANXA1 | FPR1 | 0.0258361 | 58.66667 | 1.9721433 | 16 | 0.3275053 | 28 | 0 | 132 |
| Monocytes | Monocytes | APP | APLP2 | 0.0295788 | 59.00000 | 1.9764900 | 15 | 0.3151412 | 30 | 0 | 132 |
| Monocytes | Monocytes | APP | LRP1 | 0.0295788 | 57.66667 | 1.5450213 | 30 | 0.5254924 | 11 | 0 | 132 |
| Monocytes | NK cells | S100A8 | CD69 | 0.0349511 | 64.00000 | 3.6315054 | 2 | 0.2444388 | 58 | 0 | 132 |
| B-cells | CD8+ T-cells | HLA-DRA | LAG3 | 0.0356451 | 62.66667 | 1.4660470 | 33 | 0.3762268 | 23 | 0 | 132 |
| Monocytes | Monocytes | TIMP2 | CD44 | 0.0377868 | 62.00000 | 1.8390221 | 20 | 0.2981862 | 34 | 0 | 132 |
| B-cells | CD8+ T-cells | HLA-DRB1 | LAG3 | 0.0377868 | 62.66667 | 1.4643086 | 34 | 0.3865241 | 22 | 0 | 132 |
| B-cells | B-cells | ST6GAL1 | CD22 | 0.0399878 | 63.66667 | 1.6805757 | 24 | 0.2961188 | 35 | 0 | 132 |
| Monocytes | B-cells | S100A8 | CD69 | 0.0511082 | 68.66667 | 3.6156569 | 3 | 0.2242757 | 71 | 0 | 132 |
| source | target | ligand.complex | receptor.complex | magnitude_rank | specificity_rank | magnitude_mean_rank | natmi.prod_weight | cellphonedb.lr.mean | specificity_mean_rank | logfc.logfc_comb | natmi.edge_specificity | cellphonedb.pvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Monocytes | CD8+ T-cells | S100A8 | CD69 | 0.0000076 | 0.0175097 | 1.0 | 7.326477 | 3.197897 | 62.66667 | 3.7367472 | 0.2472976 | 0 |
| Monocytes | NK cells | S100A8 | CD69 | 0.0000304 | 0.0349511 | 2.0 | 7.241781 | 3.189256 | 64.00000 | 3.6315054 | 0.2444388 | 0 |
| Monocytes | B-cells | S100A8 | CD69 | 0.0000684 | 0.0511082 | 3.0 | 6.644426 | 3.128312 | 68.66667 | 3.6156569 | 0.2242757 | 0 |
| Monocytes | CD4+ T-cells | S100A8 | CD69 | 0.0001216 | 0.0511082 | 4.0 | 5.234948 | 2.984513 | 81.33333 | 3.4042166 | 0.1767002 | 0 |
| Monocytes | Monocytes | GRN | TNFRSF1B | 0.0004864 | 0.0120408 | 6.5 | 4.482562 | 2.152250 | 52.00000 | 2.8381942 | 0.4409884 | 0 |
| Monocytes | Monocytes | ANXA1 | FPR1 | 0.0007600 | 0.0258361 | 8.5 | 3.813116 | 2.062626 | 58.66667 | 1.9721433 | 0.3275053 | 0 |
| CD8+ T-cells | CD8+ T-cells | CD48 | CD2 | 0.0012843 | 0.2995712 | 9.5 | 3.842218 | 2.015224 | 196.66667 | 0.2451358 | 0.0978671 | 0 |
| Monocytes | Monocytes | PSAP | LRP1 | 0.0012843 | 0.0056819 | 9.5 | 3.264083 | 2.194912 | 50.66667 | 2.7616439 | 0.5199923 | 0 |
| Monocytes | Monocytes | TIMP2 | CD44 | 0.0017099 | 0.0377868 | 11.5 | 3.742651 | 2.007968 | 62.00000 | 1.8390221 | 0.2981862 | 0 |
| CD4+ T-cells | CD8+ T-cells | CD48 | CD2 | 0.0019455 | 0.2703664 | 12.5 | 3.682539 | 1.963626 | 183.33333 | 0.3500985 | 0.0937999 | 0 |
| NK cells | Monocytes | ANXA1 | FPR1 | 0.0027435 | 0.0511082 | 14.5 | 3.440565 | 1.929412 | 70.66667 | 1.3956981 | 0.2955073 | 0 |
| NK cells | NK cells | CLEC2B | KLRF1 | 0.0033515 | 0.0511082 | 19.5 | 2.822475 | 1.873998 | 64.00000 | 2.0815390 | 0.2628499 | 0 |
| B-cells | CD8+ T-cells | CD48 | CD2 | 0.0036782 | 0.6492640 | 16.5 | 3.403725 | 1.873532 | 231.00000 | -0.0089090 | 0.0866980 | 0 |
| Monocytes | CD8+ T-cells | CD48 | CD2 | 0.0040202 | 0.3816175 | 17.5 | 3.339887 | 1.852904 | 208.00000 | 0.2517856 | 0.0850720 | 0 |
| CD8+ T-cells | CD4+ T-cells | CD48 | CD2 | 0.0043774 | 0.8811144 | 20.0 | 2.920353 | 1.829596 | 249.00000 | -0.0987362 | 0.0743858 | 0 |
| CD8+ T-cells | Monocytes | ANXA1 | FPR1 | 0.0047498 | 0.0511082 | 20.0 | 3.158117 | 1.828416 | 71.33333 | 1.4093857 | 0.2712481 | 0 |
| NK cells | CD8+ T-cells | CD48 | CD2 | 0.0051374 | 0.6429605 | 20.0 | 3.249439 | 1.823677 | 234.33333 | -0.0045504 | 0.0827681 | 0 |
| CD8+ T-cells | NK cells | CLEC2B | KLRF1 | 0.0055402 | 0.0511082 | 25.5 | 2.529306 | 1.819794 | 69.33333 | 2.0957862 | 0.2355478 | 0 |
| Monocytes | Monocytes | LGALS9 | CD44 | 0.0059581 | 0.0534667 | 24.0 | 2.780072 | 1.818917 | 104.00000 | 1.3670928 | 0.1556769 | 0 |
| CD4+ T-cells | CD4+ T-cells | CD48 | CD2 | 0.0068397 | 0.6059944 | 24.5 | 2.798986 | 1.777999 | 238.66667 | 0.0062265 | 0.0712944 | 0 |
| Monocytes | Monocytes | ADAM10 | CD44 | 0.0073033 | 0.1295945 | 28.0 | 2.518577 | 1.767559 | 125.66667 | 1.1033269 | 0.1231780 | 0 |
| Monocytes | NK cells | CLEC2B | KLRF1 | 0.0077821 | 0.0511082 | 30.0 | 2.198078 | 1.758553 | 76.33333 | 2.0023648 | 0.2047014 | 0 |
| B-cells | CD4+ T-cells | CD48 | CD2 | 0.0087852 | 1.0000000 | 27.5 | 2.587068 | 1.687904 | 283.33333 | -0.3527810 | 0.0658965 | 0 |
| Monocytes | B-cells | TNFSF13B | TNFRSF13C | 0.0093096 | 0.0012215 | 26.0 | 2.838560 | 1.684821 | 47.00000 | 2.7856216 | 0.8761565 | 0 |
| Monocytes | CD4+ T-cells | CD48 | CD2 | 0.0098492 | 0.8354175 | 29.5 | 2.538547 | 1.667276 | 256.00000 | -0.0920863 | 0.0646606 | 0 |
heatmap_number_celltypes = max(5, length(unique(sc$annotation)) * 0.3 + 0.5)
fig_number_celltypes = length(unique(sc$annotation)) * 3We plot the frequencies of interactions for each pair of potentially communicating cell types.
First, we filter interactions by aggregate_rank (<= 0.01), as a indicator for the robustly, highly ranked interactions concordant between methods and visualize the number of the inferred interactions between sender and retriever cell types. Here we assume that the number of the interactions is informative of the occurring communication events. However, this is a rather strong assumption affected by the applied filter criteria. Thus, any conclusions need to be confirmed by supporting information such as biological prior knowledge (Dimitrov et al. (2021)).
Heatmap displaying the interaction frequencies of all cell types.
# First, we filter interactions by aggregate_rank, as a indicator for the robustly, highly ranked interactions concordant between methods.
# These p-values are already corrected.
liana_agg_filtered = dplyr::filter(liana_agg, aggregate_rank <= param$liana_agg_rank_threshold)
number_ccinteractions = length(rownames(unique(liana_agg_filtered[,1:2])))
number_total_interactions = length(rownames(liana_agg_filtered))
# Test for sufficient number of interactions
if (number_ccinteractions*5 > number_total_interactions) {
message(paste0("The liana_agg_rank_threshold ", param$liana_agg_rank_threshold, " is too stringend. Increasing threshold automatically by factor 10."))
param$liana_agg_rank_threshold = param$liana_agg_rank_threshold * 10
# Rerun with lower filter threshold
liana_agg_filtered = dplyr::filter(liana_agg, aggregate_rank <= param$liana_agg_rank_threshold)
}×
(Message)
The liana_agg_rank_threshold 0.01 is
too stringend. Increasing threshold automatically by factor 10.
# Plot frequency heatmap
p = heat_freq(liana_agg_filtered, font_size = 10, pallette = c("white", "navy"))
pChord diagram displaying the interaction frequencies of all cell types and chosen source and target cell type groups.
p1 = wrap_elements(full= ~ chord_freq(liana_agg_filtered, cex = 0.8, facing = "inside", adj = c(0.5, 0))) +
ggtitle("All cell types")
p2 = wrap_elements(full= ~ chord_freq(liana_agg_filtered, source_groups = param$sender, target_groups = param$receiver, cex = 0.8, facing = "inside", adj = c(0.5, 0))) +
ggtitle("Selected cell types")
p = p1 + p2
pThe dotplot displays the specific ligand-receptor interactions between sender (Monocytes, NK cells) to receiver (NK cells, CD4+ T-cells, CD8+ T-cells, B-cells) cells. “Specificity” refers to how specific this interaction is to a cell type pair relative to the rest of the cell type pairs, while “Magnitude” is a direct measure of the expression alone. 30 interactions with the highest “Magnitude” and “Specificity” are shown.
p = liana_dotplot(liana_rank_agg, source_groups = param$sender, target_groups = param$receiver, ntop = 30, specificity = 'specificity_rank', magnitude = 'magnitude_rank', invert_specificity = TRUE, invert_magnitude = TRUE) +
theme(text = element_text(size = 10)) +
theme(title = element_text(size = 12, face = "bold")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1.2, colour = "black", face = "plain", size = 10)) +
theme(axis.text.y = element_text(size = 10)) +
theme(axis.title.x.bottom = element_text(vjust = -1.5)) +
theme(axis.title.y = element_text(vjust = 3.5))
p
Here, for each cell type, we plot the 10 interactions (signals) with the
highest “Magnitude” and “Specificity” that the cells receive from each
of the other cell types.
p_list = list()
for (i in unique(sc$annotation)) {
p_list[[i]] = liana_dotplot(liana_rank_agg, target_groups = i, ntop = 10, specificity = 'specificity_rank', magnitude = 'magnitude_rank', invert_specificity = TRUE, invert_magnitude = TRUE) +
AddStyle(title="", xlab=i, legend_position = "none") +
theme(text = element_text(size = 10)) +
theme(title = element_text(size = 12, face = "bold")) +
theme(axis.text.x.bottom = element_blank()) +
theme(axis.title.x = element_text(vjust = -1.5)) +
theme(axis.title.y = element_blank())
}
p = patchwork::wrap_plots(p_list, ncol = 1)
p# Add experiment details
Seurat::Misc(sc, "experiment") = list(project_id=param$project_id, date=Sys.Date(), species=gsub("_gene_ensembl", "", param$mart_dataset))
# Add parameter
Seurat::Misc(sc, "parameters") = param
# Add technical output (note: cannot use Misc function here)
sc@misc$technical = data.frame(ScrnaseqSessionInfo(param$path_to_git))### Export seurat object
saveRDS(sc, file=file.path(param$path_out, "data", "sc.rds"))### Convert data
################################################################################
### Export to Loupe Cell Browser
if (all(param$path_data$type == "10x")) {
# Export reductions (umap, pca, others)
for(r in Seurat::Reductions(sc)) {
write.table(Seurat::Embeddings(sc, reduction=r)[,1:2] %>% as.data.frame() %>% tibble::rownames_to_column(var="Barcode"),
file=file.path(param$path_out, "data", paste0("Loupe_projection_", r, ".csv")), col.names=TRUE, row.names=FALSE, quote=FALSE, sep=",")
}
# Export categorical metadata
loupe_meta = sc@meta.data
idx_keep = sapply(1:ncol(loupe_meta), function(x) !is.numeric(loupe_meta[,x]))
write.table(x=loupe_meta[, idx_keep] %>% tibble::rownames_to_column(var="barcode"),
file=file.path(param$path_out, "data", "Loupe_metadata.csv"), col.names=TRUE, row.names=FALSE, quote=FALSE, sep=",")
}
### Export as andata object (h5ad file)
# Convert Seurat v5 single cell object to anndata object
# Does not work for SCT at the moment
# However, it would contain counts and data layer
#scCustomize::as.anndata(x = sc, file_path = file.path(param$path_out, "data"), file_name = "sc_anndata.h5ad", assay="RNA",
# main_layer = "data", other_layers = "counts", transer_dimreduc = TRUE, verbose = FALSE)
# scesy only worked for v3 and v4 objects
# Convert to V3/4/Assay structure
sc_v3 = scCustomize::Convert_Assay(seurat_object = sc, convert_to = "V3", assay = "RNA")
# Convert Seurat v3 single cell object to anndata object
adata = sceasy::convertFormat(sc_v3, from="seurat", to="anndata", outFile=NULL, assay=Seurat::DefaultAssay(sc_v3))
# Write to h5ad file
adata$write(file.path(param$path_out, "data", "sc_anndata.h5ad"), compression="gzip")
### Export count matrix
invisible(ExportSeuratAssayData(sc,
dir=file.path(param$path_out, "data", "counts"),
assays=param$assay_raw,
slot="counts",
include_cell_metadata_cols=colnames(sc[[]]),
metadata_prefix=paste0(param$project_id, ".")))
### Export metadata
openxlsx::write.xlsx(x=sc[[]] %>% tibble::rownames_to_column(var="Barcode"), file=file.path(param$path_out, "data", "cell_metadata.xlsx"), rowNames=FALSE, colNames=TRUE)We export the assay data, cell metadata, clustering and visualization.
Result files:
* sc.rds: Seurat object for import into R
* cell_metadata.xlsx: Cell metadata in excel format
* counts folder: Contains raw count matrix of the entire dataset (sparse
matrix)
We export the assay data, cell metadata, clustering and visualization in andata format.
Result file:
* sc.h5ad: H5AD object
If all provided datasets are of type “10x”, we export the UMAP 2D visualization, metadata such as the cell clusters, and lists of differentially expressed genes, so you can open and work with these in the Loupe Cell Browser (https://support.10xgenomics.com/single-cell-gene-expression/software/visualization/latest/what-is-loupe-cell-browser).
Result files are:
* Loupe_projection_(umap|pca|…).csv: Seurat UMAP/PCA/… projections for
visualization
* Loupe_metadata.csv: Seurat cell meta data including clusters and cell
cycle phases
Projections can be imported in Loupe via “Import Projection”, cell meta data via “Import Categories” and gene sets via “Import Lists”.
The following parameters were used to run the workflow.
out = ScrnaseqParamsInfo(params=param)
knitr::kable(out, align="l") %>%
kableExtra::kable_styling(bootstrap_options=c("striped", "hover"), full_width=FALSE, position="left")| Name | Value |
|---|---|
| path_to_git | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis |
| path_to_basic_settings | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/basic_settings.R |
| path_to_advanced_settings | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/advanced_settings.R |
| scriptname | scripts/ccc_analysis/ccc_analysis.Rmd |
| author | kosankem |
| assay_raw | RNA |
| downsample_cells_equally | FALSE |
| cell_filter | nFeature_RNA:20, NA; nCount_RNA:200, 20000; percent_mt:0, 18; Sample1:nFeature_RNA=c(20, NA), nCount_RNA=c(200, 20000), percent_mt=c(0, 18); Sample2:nFeature_RNA=c(20, NA), nCount_RNA=c(200, 20000), percent_mt=c(0, 18) |
| feature_filter | min_counts:1; min_cells:5; Sample1:min_counts=1, min_cells=5; Sample2:min_counts=1, min_cells=5 |
| samples_min_cells | 10 |
| norm | RNA |
| cc_remove | FALSE |
| cc_remove_all | FALSE |
| cc_rescore_after_merge | TRUE |
| integrate_samples | method:merge; dimensions:30; k_anchor:10; k_weight:100; reference:; integration_function:CCAIntegration |
| experimental_groups | homogene |
| pc_n | 8 |
| cluster_k | 20 |
| umap_k | 30 |
| cluster_resolution_test | 0.5, 0.7, 0.8 |
| cluster_resolution | 0.6 |
| marker_padj | 0.05 |
| marker_log2FC | 1 |
| marker_pct | 0.25 |
| enrichr_site | Enrichr |
| enrichr_padj | 0.05 |
| col_palette_samples | ggsci::pal_igv |
| col_palette_clusters | ggsci::pal_igv |
| col_palette_annotation | ggsci::pal_ucscgb |
| col | #0086b3 |
| col_bg | #D3D3D3 |
| pt_size | 0.5 |
| liana_methods | logfc, natmi, cellphonedb |
| liana_agg_rank_threshold | 0.1 |
| project_id | Testdata |
| data | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/output/Testdata/cluster_analysis/data/sc.rds |
| path_data | name:Sample1, Sample2; type:10x, 10x; path:/mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/data/counts/sample1/, /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/data/counts/sample2/ |
| species | human |
| sender | Monocytes, NK cells |
| receiver | NK cells, CD4+ T-cells, CD8+ T-cells, B-cells |
| path_out | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/output/Testdata/ccc_analysis |
| debugging_mode | default_debugging |
| mart_dataset | hsapiens_gene_ensembl |
| mt | ^MT- |
| enrichr_dbs | GO_Biological_Process_2023, WikiPathway_2023_Human, Azimuth_2023, CellMarker_2024 |
| annotation_dbs | BlueprintEncodeData() |
| annot_version | 98 |
| annot_main | ensembl=ensembl_gene_id, symbol=external_gene_name, entrez=entrezgene_accession |
| mart_attributes | ensembl_gene_id, external_gene_name, entrezgene_accession, chromosome_name, start_position, end_position, percentage_gene_gc_content, gene_biotype, strand, description |
| path_reference | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/references/hsapiens_gene_ensembl/98 |
| reference | hsapiens_gene_ensembl.v98.annot.txt |
| file_annot | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/references/hsapiens_gene_ensembl/98/hsapiens_gene_ensembl.v98.annot.txt |
| file_cc_genes | /mnt/ngsnfs/single_cell_dev/scRNAseq_processing/sc_analysis/references/hsapiens_gene_ensembl/98/cell_cycle_markers.xlsx |
| col_samples | Sample1=#5050FFFF, Sample2=#CE3D32FF |
| col_clusters | 1=#5050FFFF, 2=#CE3D32FF, 3=#749B58FF, 4=#F0E685FF, 5=#466983FF, 6=#BA6338FF, 7=#5DB1DDFF |
| col_annotation | Monocytes=#FF0000FF, CD8+ T-cells=#FF9900FF, CD4+ T-cells=#FFCC00FF, B-cells=#00FF00FF, NK cells=#6699FFFF |
This report was generated using the scripts/ccc_analysis/ccc_analysis.Rmd script. Software versions were collected at run time.
out = ScrnaseqSessionInfo(param$path_to_git)
knitr::kable(out, align="l") %>%
kableExtra::kable_styling(bootstrap_options=c("striped", "hover"))| Name | Value |
|---|---|
| Run on: | Wed Sep 04 11:16:49 AM 2024 |
| ktrns/scrnaseq | 8a210a92c96c41377c5123e641ec837d56bf06ae |
| Container | NA |
| R | R version 4.2.1 (2022-06-23) |
| Platform | x86_64-pc-linux-gnu (64-bit) |
| Operating system | Ubuntu 20.04.6 LTS |
| Host name | hpc-rc11 |
| Host OS | #138-Ubuntu SMP Wed Jun 22 15:00:31 UTC 2022 (5.4.0-122-generic) |
| Packages | abind1.4-5, backports1.4.1, basilisk1.17.2, basilisk.utils1.17.2, beachmat2.14.2, beeswarm0.4.0, bibtex0.5.1, Biobase2.58.0, BiocGenerics0.44.0, BiocNeighbors1.16.0, BiocParallel1.32.6, BiocSingular1.14.0, bit4.0.5, bit644.0.5, bitops1.0-7, blob1.2.4, bluster1.8.0, bslib0.6.1, cachem1.0.8, Cairo1.6-2, cellranger1.1.0, checkmate2.3.1, circlize0.4.16, cli3.6.2, clue0.3-65, cluster2.1.3, clustree0.5.1, codetools0.2-18, colorspace2.1-0, ComplexHeatmap2.14.0, cowplot1.1.3, crayon1.5.2, curl5.2.1, data.table1.15.2, DBI1.2.2, DelayedArray0.24.0, DelayedMatrixStats1.20.0, deldir2.0-4, digest0.6.34, dir.expiry1.6.0, doParallel1.0.17, dotCall641.1-1, dplyr1.1.4, dqrng0.3.2, edgeR3.40.2, ellipsis0.3.2, enrichR3.2, evaluate0.23, fansi1.0.6, farver2.1.1, fastDummies1.7.3, fastmap1.1.1, filelock1.0.3, fitdistrplus1.1-11, forcats1.0.0, foreach1.5.2, future1.33.1, future.apply1.11.1, generics0.1.3, GenomeInfoDb1.34.9, GenomeInfoDbData1.2.9, GenomicRanges1.50.2, GetoptLong1.0.5, ggbeeswarm0.7.2, ggforce0.4.2, ggplot23.5.0, ggprism1.0.5, ggraph2.2.1, ggrastr1.0.2, ggrepel0.9.5, ggridges0.5.6, ggsci3.0.1, GlobalOptions0.1.2, globals0.16.2, glue1.7.0, goftest1.2-3, graphlayouts1.1.1, gridExtra2.3, gridGraphics0.5-1, gtable0.3.4, highr0.10, hms1.1.3, htmltools0.5.7, htmlwidgets1.6.4, httpuv1.6.14, httr1.4.7, ica1.0-3, igraph2.0.2, IRanges2.32.0, irlba2.3.5.1, iterators1.0.14, janitor2.2.0, jquerylib0.1.4, jsonlite1.8.8, kableExtra1.4.0, KernSmooth2.23-20, knitcitations1.0.12, knitr1.45, labeling0.4.3, later1.3.2, lattice0.20-45, lazyeval0.2.2, leiden0.4.3.1, liana0.1.14, lifecycle1.0.4, limma3.54.2, listenv0.9.1, lmtest0.9-40, locfit1.5-9.9, logger0.3.0, lubridate1.9.3, magrittr2.0.3, MASS7.3-58, Matrix1.6-5, MatrixGenerics1.10.0, matrixStats1.1.0, memoise2.0.1, metapod1.6.0, mime0.12, miniUI0.1.1.1, munsell0.5.0, nlme3.1-157, OmnipathR3.13.22, openxlsx4.2.5.2, paletteer1.6.0, parallelly1.37.1, patchwork1.2.0, pbapply1.7-2, pillar1.9.0, pkgconfig2.0.3, plotly4.10.4, plyr1.8.9, png0.1-8, polyclip1.10-6, prettyunits1.2.0, progress1.2.3, progressr0.14.0, promises1.2.1, purrr1.0.2, R.methodsS31.8.2, R.oo1.26.0, R.utils2.12.3, R62.5.1, ragg1.2.7, RANN2.6.1, rappdirs0.3.3, RColorBrewer1.1-3, Rcpp1.0.12, RcppAnnoy0.0.22, RcppHNSW0.6.0, RCurl1.98-1.14, readr2.1.5, readxl1.4.3, RefManageR1.4.0, rematch22.1.2, remotes2.4.2.1, renv0.16.0, reshape21.4.4, reticulate1.35.0, rjson0.2.21, rlang1.1.3, rmarkdown2.25, ROCR1.0-11, RSpectra0.16-1, RSQLite2.3.5, rstudioapi0.15.0, rsvd1.0.5, Rtsne0.17, rvest1.0.4, S4Vectors0.36.2, sass0.4.8, ScaledMatrix1.6.0, scales1.3.0, scater1.26.1, scattermore1.2, scCustomize2.1.2, sceasy0.0.7, scran1.26.2, sctransform0.4.1, scuttle1.8.4, sessioninfo1.2.2, Seurat5.0.2, SeuratObject5.0.1, shape1.4.6.1, shiny1.8.0, SingleCellExperiment1.20.1, snakecase0.11.1, sp2.1-3, spam2.10-0, sparseMatrixStats1.10.0, spatstat.data3.0-4, spatstat.explore3.2-6, spatstat.geom3.2-9, spatstat.random3.2-3, spatstat.sparse3.0-3, spatstat.utils3.0-4, statmod1.5.0, stringi1.8.3, stringr1.5.1, SummarizedExperiment1.28.0, survival3.2-13, svglite2.1.3, systemfonts1.0.5, tensor1.5, textshaping0.3.7, tibble3.2.1, tidygraph1.3.1, tidyr1.3.1, tidyselect1.2.0, timechange0.3.0, tweenr2.0.3, tzdb0.4.0, utf81.2.4, uwot0.1.16, vctrs0.6.5, vipor0.4.7, viridis0.6.5, viridisLite0.4.2, withr3.0.0, WriteXLS6.7.0, xfun0.41, XML3.99-0.16.1, xml21.3.6, xtable1.8-4, XVector0.38.0, yaml2.3.8, zip2.3.1, zlibbioc1.44.0, zoo1.8-12 |
This Workflow was developed as module of the sc_analysis workflow at the Research Core Unit Genomics, Hannover Medical School (Hannover, Germany). Seurat Vignettes were initially used as templates (Hao et al. (2023), Hao et al. (2021)).
# Writes knitcitations references to references.bib file.
knitcitations::write.bibtex(file=file.path(param$path_out, "references.bib"))